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It is shown that the nonequilibrium self-energy of an interacting lattice-fermion model has a unique Lehmann 
representation. Based on the construction of a suitable non-interacting effective medium, we provide an ex¬ 
plicit and numerically practicable scheme to construct the Lehmann representation for the self-energy, given the 
Lehmann representation of the single-particle nonequilibrium Green’s function. This is of particular importance 
for an efficient numerical solution of Dyson’s equation in the context of approximations where the self-energy is 
obtained from a reference system with a small Hilbert space. As compared to conventional techniques to solve 
Dyson’s equation on the Keldysh contour, the effective-medium approach allows to reach a maximum propa¬ 
gation time which can be several orders of magnitude longer. This is demonstrated explicitly by choosing the 
nonequilibrium cluster-perturbation theory as a simple approach to study the long-time dynamics of an inhomo¬ 
geneous initial state after a quantum quench in the Hubbard model on a 10 x 10 square lattice. We demonstrate 
that the violation of conservation laws is moderate for weak Hubbard interaction and that the cluster approach 
is able to describe prethermalization physics. 

PACS numbers: 71.10.-w,71.10.Fd,67.85.Lm,78.47.J- 


I. INTRODUCTION 

The study of physical phenomena that arise in strongly 
correlated systems far from equilibrium has become a field 
of highly active research recently.DEl For the theoretical de¬ 
scription of such systems. Green’s-function-based approaches 
starting from the Keldysh formalistrPl have proven to be very 
useful. A number of different approximation schemes rely on 
this concept.EK2uentral to these approaches is the self-energy 
which is related to the one-particle Green’s function through 
Dyson’s equation. However, while the numerical solution of 
Dyson’s equation is rather straightforward in the equilibrium 
case, the computational effort is considerably increased for 
systems out of equilibrium since operations with matrices de¬ 
pending on two independent contour time variables typically 
scale cubically in the number of time steps. Apart from other 
challenges characteristic for the respective approach, already 
this scaling poses a severe limit on the maximal reachable 
propagation time in a numerical calculation. Applying ad¬ 
ditional concepts or a pprox imations, such as the generalized 
Kadanoff-Baym ansat J^IHI or exploiting a rapid decay of the 
memory,G31 are necessary to overcome this limitation. 

It was proposed recentlyC^! that it can be advantageous to 
avoid the direct inversion of Dyson’s equation by applying a 
mapping onto a Markovian propagation scheme. To this end 
it is necessary to assume the existence of a certain functional 
form for the nonequilibrium self-energy, namely the existence 
of a Lehmann representation. 

In the present paper we explicitly construct this Lehmann 
representation. With this at hand, we pick up the proposed 
idea to solve Dyson’s equation by means of a Markovian prop¬ 
agation and exploit the fact that the Lehmann representation 
of the exact self-energy of a small reference system has a fi¬ 
nite number of terms only. This allows us to solve Dyson’s 
equation with an effort that scales linearly in the maximum 
propagation time fmax- 

For equilibrium Green’s functions^e Lehmann represen¬ 
tation is a well established concept.^I^ It uncovers the ana¬ 


lytical properties of the Green’s function and can be used 
to show that the related spectral function is positive defi¬ 
nite. It is further essential for the evaluation of diagrams 
through contour integrations in the complex frequency plane, 
for the derivation of sum rules, etc. The generalization of 
the Lehmann representation to nonequilibrium Green’s func¬ 
tion is straightforward.El Applications include nonequilibrium 
dynamical mean-field theory (DMFT) where it allows for a 
Hamiltonian-based formulation of the impurity problem.El 

The explicit construction of a Lehmann representation for 
the self-energy, on the other hand, turns out to be more te¬ 
dious, already for the equilibrium case; In a recent work such 
a construction was worked ouP^ from a diagrammatic per¬ 
spective and used to cure the problem of possibly negative 
spectral functions arising from a summation of a subclass of 
diagrams. 

Here, we address the nonequilibrium self-energy of a gen¬ 
eral, interacting lattice-fermion model; (i) We rigorously 
show the existence of the Lehmann representation by pre¬ 
senting an explicit construction scheme that is based on 
the Lehmann representation of the nonequilibrium Green’s 
function, (ii) Usi ng a simp le example, namely the cluster- 
perturbation theor}Ll21i2HIir (CPT), we furthermore demon¬ 
strate that the Lehmann representation of the self-energy can 
in fact be implemented numerically and used to study the time 
evolution of a locally perturbed Hubbard model on a large 
square lattice (10 x 10 sites). Propagation times of several 
orders of magnitude in units of the inverse hopping ampli¬ 
tude can be reached with modest computational resources, 
(iii) While the CPT approximation for the self-energy is rather 
crude and shown to violate a number of conservation laws, it is 
possible with this approximation to study the weak-coupling 
limit of the Hubbard model in a reasonable way. In particu¬ 
lar we demonstrate that prethermalization physics is already 
captured on this level. 

The paper is organized as follows; In Section|I^we briefly 
discuss the generalization of the Lehmann representation to 
nonequilibrium Green’s functions. The main idea of Ref. ITSl 
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about the Markovian propagation scheme is recalled in Sec. 
|111 A| The explicit construction scheme for the nonequilibrium 
self-energy is outlined in Section 111 B Section IV is devoted 
to the application of our formalism to the cluster-perturbation 
theory. Sec. [^presents numerical results for the time evolu¬ 
tion of a local perturbation in the fermionic Hubbard model. 
We conclude the paper with a summary and an outlook in Sec. 

m 


II. LEHMANN REPRESENTATION OF THE 
ONE-PARTICLE GREEN’S FUNCTION 
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FIG. 1: Keldysh-Matsubara contour C. Ci denotes the upper branch, 
C 2 the lower branch and C 3 the Matsubara branch. In the shown 
example t is later than t' in sense of the contour, denoted as t >c t' 
in the text. 


We consider an arbitrary, fermionic model Hamiltonian 


H{t) = Y,{Tij{t) - 8ijl^)c]cj -b 
ij 


^ iji'f 


( 1 ) 


where the indices i,j run over the possible one-particle 
orbitals (lattice sites, local orbitals, spin projection, ...). 
Fermions in such states are created (annihilated) by the op¬ 
erators cj (c,). At time t — 0, the system with Hamilto¬ 
nian H{Q) = Mini is assumed to be in thermal equilibrium 
with inverse temperature j3 and chemical potential /r. Non¬ 
equilibrium real-time dynamics for f > 0 is initiated by the 
time dependence of the one-particle or the interaction param¬ 
eters. This covers challenging experimental setups such as 
time-resolved photoemission spectroscopjEU or experiments 
with ultracold gases in optical lattices.lSl 

The one-particle Green’s function is given by 


= -i{TcCi{t)c]{t'))H 



(exp{-j5Hini) Tcci{t)c]{t') y 


( 2 ) 


Here, /(e) = + 1) ^ denotes the Fermi-function while 

refers to the contour variant of the Heaviside step 
function (0c(f,fO = 1 for f >c t', = 0 otherwise). 

Q{t) is defined to be equal on the upper and lower branch of 
the contour and furthermore constant on the Matsubara branch 
with Q{—iT) = 2(0) and T S [0,j3]. If the eigenstates \m) of 
the initial Hamiltonian (i.e., Hini\ni) = Em\m)) are used as a 
basis for tracing over the Fock space in Eq. (|^, one has 

Qia(t) = Qi{m,n){t) = , (5) 

where zp„.n) = V+ e-P^")/Z and where the su¬ 
perindex a = {m,n) labels the possible one-particle 
excitations with corresponding excitation energies 
Ea = £(m,n) = E„ — Em. Note that this definition of Q{t) 
indeed satisfies Q{—iz) = 2(0). We emphasize that 2 as a 
matrix is not quadratic. Our expression can be seen as a direct 
generalization of the time-independent 2-matrix discussed in 
Ref |27] We further note that the rows of the 2-matrix fulfill 
the orthonormality condition 


where “tr(...)” traces over the Fock space, i.e., we 
take averages using the grand-canonical ensemble. 
Z = tr(exp(—j3//ini)) defines the grand-canonical parti¬ 
tion function and 7c the time-ordering operator on the 
L-shaped Keldysh-Matsubara contour C (see Fig. [^l. The 
time variables t and t' are understood as contour times 
that can lie on the upper, lower or Matsubara branch 
of C. We further introduce the convention that opera¬ 
tors with a hat carry a time dependence according to the 
Heisenberg picture, i.e., c,(f) = t/^(f,0)c,t/(f,0), where 
U{t,t') = Texp (—//(//(fi)dfi) is the system’s time- 
evolution operator and T the time-ordering operator. An 
in-depth introduction to the Keldysh formalisrrPlcan be found 
in Refs.|25l|26l 

As has been shown in Ref. [m the one-particle Green’s 
function can be cast into the form 

=Y,Qia{t)gi£a;t/)Q%it’), (3) 

a 

which we will call its Lehmann representation in the follow¬ 
ing. g(e;f,f') is the non-interacting Green’s function of an 
isolated one-particle mode (/zmode = ec^c) with excitation en¬ 
ergy e: 


m)QHt)]ij =E2,a(f)ej„(o = (|cMo,c](o})// = 5 ,;, 

(6) 

where {A,B} =AB + BA denotes the anticommutator. 


HI. LEHMANN REPRESENTATION OF THE 
SELF-ENERGY 

A. Motivation 

In several Green’s-function-based methods, an approximate 
self-energy E' is obtained from a small reference system us¬ 
ing exact diagonalization. The desired one-particle Green’s 
function G of a much larger system is then obtained through 
Dyson’s equation 

Gijit/) = [Go]ij{t/) (7) 

where Go denotes the non-interacting Green’s function (i.e., 
G = 0) of the model given by Eq. 0. Typical exam¬ 
ples include dynamical mean-field theory (DMET),E02lIlll 
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where E' is obtained from a single -impurity Anderson 
modeljSlor cluster-p erturbatiorP^^^Jt^ and self-energy func¬ 
tional theory pSl^llll where E' stems from a small reference 
system. To solve Eq. 0 numerically, a discretization of the 
continuous time-contour C is necessary. The number of time 
steps required to reach a given maximal time is dependent on 
the lowest relevant timescale that is set by a given Hamilto¬ 
nian. Based on this discretization, the effort required to solve 
Eq. (0 for G scales cubically in the number of time steps and 
also the system size. Despite this challenge also the mem¬ 
ory consumption, which scales quadratically in these quanti¬ 
ties, poses a problem. Progress was made recentljG^by intro¬ 
ducing a mapping of Eq. Q onto a Markovian propagation- 
scheme. 

The idea proposed by the authors of Ref. [15] relies on the 
assumption that the self-energy can be written in the following 
form: 

E;/f,f') = (8) 

Here, E'™(f) denotes the time-local Hartree-Eock term. This 
decomposition is very similar to the expression Eq. ([^ for the 
Green’s function. We will refer to this as the Lehmann repre¬ 
sentation of the self-energy. The immediate and important ad¬ 
vantage of the Lehmann representation is tha t the s elf-energy 
can be interpreted as a hybridization function.^ElThis prop¬ 
erty allows to write down an effective non-interacting model 
with Hamiltonian 

u 

+ Y,ihis{t)c]as + h.c.) +Y^hssalas. 

is s 

The s-degrees of freedom represent “virtual” orbitals in ad¬ 
dition to the physical degrees of freedom labeled by i. They 
form an “effective medium” with on-site energies hss and hy¬ 
bridization strengths such that the interacting Green’s 

function of the original model is the same as the Green’s func¬ 
tion of the effective non-interacting model on the physical or¬ 
bitals: 

Gijit/) = (10) 

With this simple construction, the inversion of the Dyson 
equation can be avoided in favor of a Markovian time prop¬ 
agation within a non-interacting model. 

As a successful benchmark, an interaction quench in an in¬ 
homogeneous Hubbard model was treated with nonequilib¬ 
rium DMET in Ref. [TSjusing self-consistent second-order per¬ 
turbation theory as impurity solver. On the theoretical side, 
however, it remained an open question if the existence of a 
Lehmann representation must be postulated or if this is a gen¬ 
eral property of the nonequilibrium self-energy. 

In the following we explicitly derive Eq. 0 for the exact 
self-energy corresponding to the general, interacting Hamilto¬ 
nian defined in Eq. 0,i.e., we show that the exact self-energy 
can always be written in the form of a Lehmann representa¬ 
tion. The proposed construction scheme is not only useful 
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FIG. 2: Unitary completion of the time-dependent Matrix Q{t). 
The matrix Q^{t) contains a completing set of orthonormal ba¬ 
sis vectors in its rows. For convenience, the phase factor 
£aa'{t) = 5aa'exp(—ifiat) is also absorbed into 0{t). The gener¬ 
ating, Hermitian matrix h{t) (cf. Eq. ( |12| l) can be assumed to be di¬ 
agonal in the virtual sector. 

as an analytical tool but also well suited for numerical appli¬ 
cations where an approximate self-energy is obtained from a 
small reference system using exact diagonalization. In this 
case the number of virtual orbitals is constant and the effort 
for solving Eq. Q scales linearly in fmax- This is a great ad¬ 
vantage if one is interested in long-time dynamics. 

B. Explicit construction 

We start our construction from the Lehmann representa¬ 
tion of G as stated in Eq. 0. Eor our model Hamiltonian 
o the associated one-particle excitation energies £« and the 
2-matrix are given by Eq. 0. The self-energy is related to 
this representation through Dyson’s equation E = Gq ' — G 
However, the inverse G * cannot directly be calculated with 
Eq. 0 since Q{t) is not quadratic. As a first step we block up 
the matrix Q{t) to a quadratic form. This is achieved by in¬ 
terpreting its orthonormal rows (cf. Eq. 0) as an incomplete 
set of basis vectors. Q{t) itself is an incomplete unitary trans¬ 
form from this viewpoint. We now pick an arbitrary, pairwise 
orthonormal completion of this basis to find an unitary trans¬ 
form 0{t) that contains Q{t) in its upper block (cf. Eig. 0. 
The next steps of our discussion will be independent of the 
particular completion that is chosen. The only mathematical 
requirement is that it is as smooth (and thus differentiable) 
in the time variable t as 2(f); see Appendix [ a| for numerical 
details on the construction of 0{t). 

The completed unitary transform 0{t) describes additional 
virtual orbitals (labeled by the index s, see Eig.0and Eq. 0). 
Eor convenience, we also absorb in the definition of G(f) the 
extra factor £aa'{t) = that stems from the 

non-interacting Green’s function g{ea',t,t) (cf. Eqs. 0 and 
©)■ Eor clarity in the notations we use the following index 
convention throughout this paper 

physical orbitals: i,j, virtual orbitals: r,s, 
physical or virtual orbitals: x,y, excitations: a, a'. (11) 

Like every time-dependent unitary transform, 2(f) is gener¬ 
ated by an associated Hermitian matrix. We define 

h,y{t) =Y,[idrO,ait)]Oiy{t). (12) 

a 
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Indeed, by integration we have 

6)(f) =rexp (^-i (9(0) (13) 

and furthermore h{t) is Hermitian; 

h{t) = [id,0{t)]0^(t) = idt[0{t)0'{t)] —0{t)idt0"'^(t) 

= ([idrOit)]OHt)y =hHt). (14) 

We now require the virtual part to be diagonal and 

time-independent, i.e., = hss{Q)5ss'- To this end we use 

our freedom in choosing the completing basis vectors Q^{t) 
which allows us to perform the associated unitary transform 
in the virtual sector (see Fig.l^. With the resulting hxy{t) we 
define the single-particle Harmltonian Heff{t) 

= ^^hxy(t)c^Cy, (15) 

-ty 


completion of ()(?) ^ unitary transform (9(f) leads to a valid 

effective Hamiltonian. The physical sectors of (9(f) and h(t) 
remain independent of its choice. The virtual sectors, on the 
other hand, are affected and only the special choice of (9(f) 
(cf the discussion above and below Eq. ( fTSj l) guarantees a di¬ 
agonal form of the effective medium. 

Having found an effective, non-interacting model that re¬ 
produces the correct Green’s function, it remains to link this 
back to the self-energy. The time-non-local (correlated) part 
T,fj{t,t') follows by tracing out the virtual orbitals. This proce¬ 
dure is straightforward as they are all non-interacting and we 
can use, e.g., a cavity-like ansatJl^ or an equation of motion 
based approach.C^This results in a hybridization-like function 

ET(f,f') = Y,hixit)gihxx-,t/)h%it’) (21) 

A' 

that encodes the influence of the virtual sites on the physical 
sector. The Green’s function at the physical orbitals is then 
obtained from a Dyson-like equation 


which has precisely the form of the effective Hamiltonian 
stated in Eq. (|^. The requirement of a diagonal virtual sector 
defines the effective Hamiltonian uniquely up to rotations in 
invariant subspaces. 

At time f = 0, the effective medium can be stated in a diag¬ 
onal form which is useful for the evaluation of the correspond¬ 
ing one-particle Green’s function. We recall that we required 
(9(f) to be as smooth as Q{t) and take a look at 

HO{t)]r=o = /((O)G(O) = G(0)M, (16) 

where M = (9l(0)/!(0)(9(0). Eq. ([T^ implies in particular that 
[idtQ{t)£{t)]t=Qi = Q{0)M (cf. Eig.|^. However, fromEq. (|^ 
one easily evaluates [[idtQ{t)£{t)]^t=o = Q{0)ia£a and we 
can thus identify = SaaiEa- Putting everything together 
we find 

hxyiO) = '£Oxa{0)eaO;„iO). (17) 

a 

We require that the effective medium is initially in thermal 
equilibrium with the same inverse temperature j3 and the same 
chemical potential jj. as the physical system. The associated 
one-particle Green’s function of the medium is defined as 

Fxy{t,t') = -i{TcCx{t)cl{t'))H^^f (18) 

Recalling the diagonal form of the effective medium at f = 0 
(cf. Eq. [TT] ) and using that the effective Hamiltonian ( [T5| l is 
non-interacting, we can easily rewrite this expression into 


Fij{t,t') = 


TT'-EC 


where 


{t/), 




( 22 ) 


[FQ%j{tA = [idt-hij{t)]5c{t,t'), (23) 

with 5c{t,£) = dt&c{t,£) as the contour delta function. 

To make the final connection to the self-energy we evaluate 
the physical sector of h. With 

idt Qi{m,n){t)e^* =Z(m.n){m\[ciit),H{t)]\n) 
j 

+ E (24) 

ji'f 

we obtain 


hijit) = Tijit) - -t-E|]^(f), 

Ef (f) = 2^t/,,,.,.,(f)(rc4(f)c>(f)>//e,f (25) 

i'f 


At the physical orbitals the effective Hamiltonian is thus de¬ 
termined by the Hartree-Eock Hamiltonian. By comparison of 
Eq. with the Dyson equation 


where 







(26) 


Fxy{t,t') = iY^Oxam{£a)-@c{t,t’)]0;.^{t). (19) 

a 

The physical sector of F is by construction identical with the 
Lehmann representation of G: 

i'y(f,f') =Eaa(fk(ea;fT')eJa(f') = G,7(f,f')- (20) 

a 

F encodes the full information on the one-particle excitations 
of the system defined by the Hamiltonian (0. Eq. ( |20l i fur¬ 
ther stresses the fact that in principle any (sufficiently smooth) 


[GQ%{t,t') = [id, - 5c{t,t'), (27) 

we finally identify 

E,7(f,f') = (28) 

concluding our construction of the self-energy. Let us stress 
that with Eqs. ([T^, ([^ and ( |25] l we now have an explicit 
recipe to construct the Lehmann representation of the self¬ 
energy. This representation is further unique as follows from 
the uniqueness of the corresponding effective Hamiltonian (cf. 
the discussion above and below Eq. ([T5]l). 
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C. Useful properties 

With the Hamiltonian of the effective medium, Eq. ( fTSl l, at 
hand, a number of useful properties follow immediately; 

1. Positive spectral weight 

By taking a look at the Matsubara branch only, one can 
link the Lehmann representation of the self-energy to the pos¬ 
itive dehniteness of its equilibrium spectral function. With 
E^(t — t') = —!E(—/t, — !t') we can perform the usual 
Fourier transform from imaginary time to Matsubara frequen¬ 
cies and then hnd the analytical continuation E^(a)) to the 
complex-frequency plane (see for example Ref. il7]l. The 
spectral function is dehned as 

Cfj(co) = ^ [Jlfjico + iO) - - iO)] (29) 

for real co. This can explicitly be calculated from the parame¬ 
ters of the effective Hamiltonian. One hnds: 

Cfj{(o) = '£hi,iO)h*M5ico-hss), (30) 

where 5(a)) is the Dirac delta function. The positive dehnite¬ 
ness for every O) is immediately evident. 


where £« are the excitations energies of //jni. The 6>-matrix is 
continuous at f = 0 despite the quantum quench (it only de¬ 
pends on c,(f), cf. Eq. Q). Its time derivative, however, is not 
and thus h{t) jumps from h{Q) to 

= '^[idrO,a{t)Uo^O%{Q). (32) 

a 

After this jump, the effective Hamiltonian will in general not 
be constant for times t > 0, i.e., h{t) 7 ^ /)(0+). 


IV. APPLICATION TO CLUSTER-PERTURBATION 
THEORY 

The simplest numerical applicat ion o f our formalism is 
given by cluster-perturbation theorjEMlES (CRT). The idea 
of CRT is to split the system into small clusters which can 
be treated by means of exact-diagonalization techniques. The 
cluster self-energies are then used as approximate input for 
the Dyson equation 0 to obtain the CRT Green’s func¬ 
tion. The saine co ncept is part of more powerful appro aches 
like DMFT®2S1221 or self-energy functional theotyG^BO] where 
the CRT Green’s function is self-consistently or variationally 
linked to the self-energy of a reference system. The following 
construction of an effective Hamiltonian for CRT applies to 
such techniques as well. 


2. Higher-order correlation functions 

The self-energy and its time derivatives can be used to cal¬ 
culate certain expectation values of higher order. Rrominent 
examples include the interaction energy or the local double 
occupation. Their calculation is based on the evaluation of 
contour integrals of the form /(-ck'E(f,f')G(f',f)- By com¬ 
paring the equations of motion for Gijitf') and Fxy{t,t') one 
readily hnds the identity 

j J 

+ Y^hi,{t)F,i,{t,t'). (31) 


A. Cluster-perturbation theory (CPT) 


From now on we restrict ourselves to the fermionic Hub¬ 
bard model. The locality of its interaction term allows us to 
cast its Hamiltonian into the following form; 


H{t)='£ 

1 


- lddij]cl^cija + U{t)Y,nipnjii 

ijo i 


cluster system Hj 

IfJijrr 


(33) 


inter-cluster hopping 


This is a remarkable relation as the contour integration can be 
avoided in favor of a simple matrix multiplication. 


3. Quantum quenches 

A convenient tool to drive quantum systems out of equi¬ 
librium is given by the so-called quantum quenches. Here, 
one (or more) parameters of the system are changed sud¬ 
denly. This sudden change rehects itself as a discontinuous 
time dependence of the effective Hamiltonian; Assume that 
the system is subjected to a quench at time f = 0, so that 
Hmi = const. Initially the system is in thermal equi¬ 

librium and the effective Hamiltonian is given by Eq. 


I 

J..... 


T’// 

jj 

'J'lj 

if 


J 




FIG. 3; Illustration of the partitioning of an infinite, two-dimensional 
square lattice into 2x2 clusters. The sites i, j lie within the same 
cluster /, / belongs to a different cluster J. The cluster diagonal 
part of the hopping matrix pF describes the intra-cluster, the cluster 
off-diagonal part tF, (I J) the inter-cluster-hopping. 
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Here, the indices 7,7 label the cluster systems, while the in¬ 
dices i,j run over the sites within a cluster only (see Fig. [^. 
Of course, this is fully equivalent with the usual form of the 
Hubbard model which is re-obtained by combining {I,i) to a 
superindex, i.e., {I,i) —^ i. The operator n/,<j = c\^^cjia mea¬ 
sures the particle density with spin projection a =t)i- The 
Green’s function of the isolated cluster I with intra-cluster 
Hamiltonian Hi is 

= -i{Tcdnait)c].,(t'))H,. (34) 

so that 


(35) 

where YJ denotes the corresponding self-energy. We further 
define 


G' = 


(G^ 

0 



V : 



(36) 


With the inter-cluster (ic) hopping [7’“^];^^ = (1 ~ 
CPT Green’s function is defined as 


pCPT _ ^ 

“ (G')^* -7’“ 


1 


Gn'-j:' 


(37) 


where [Gq = [id, - - iidijdij)]dc{t,t'). The 

definition of G*'^^ reveals that CPT becomes exact in the 
limit of vanishing interaction. We then have E' = 0 and thus 
GCpt = [Gq = Go. Solving Eq. ( |J7l i in case of non¬ 
vanishing E', on the other hand, requires the solution of a 
Dyson equation. This brings us back to our original problem. 


B. Application of the Lehmann representation for the 
self-energy 


Using our results from Sec. Ill we can avoid the solution of 
the Dyson equation and rather decompose the self-energies of 
the isolated clusters into their Lehmann representations: 


+ L hi, {t)8K,-,t/) [h%, (f'). (38) 


Here, h\t) are the parameters of the effective medium corre¬ 
sponding to the 7-th cluster. We define 


h\t) 


fh\t) 0 

0 h^{t) 


(39) 


It is now straightforward to realize that the inclusion of the 
inter-cluster hopping by means of Eq. is completely triv¬ 
ial in this language. Namely, 


( 40 ) 


With 


(41) 

IJ xya 

we then have 

= -i{TcCiiait)c]j,{t'))fjcPT. (42) 

While this is an easy and intuitive description, we re¬ 
mark that H{t) includes virtual orbitals. The inter-cluster 
hopping 7’“^(f), on the other hand, is defined solely in 
the physical sector and has to be blocked up accordingly 

=0). 

As an important observable we briefly discuss the calcula¬ 
tion of the total energy within CPT. While the kinetic energy 
follows straightforwardly from the one-particle density matrix 
as £^ 11(0 = -iLuLija T//aG‘-'j,{t,t+), the interaction energy 
can only be accessed indirectly through the self-energy. It is 
given by 

£mt(f) = -lE / dfiE^^(f,fi)GCPT(fi,f+). (43) 
i] 

The evaluation of this contour-integral in E q. (|43] l is straight¬ 
forward within our formalism by using Eq. (|3lt. 

V. NUMERICAL RESULTS 
A. Prethermalization 

The study of real-time dynamics initiated by an interac¬ 
tion quench in the Hubbard model has attracted much atten¬ 
tion recently.l^^HMI Here, the system is prepared in a thermal 
(usually non-interacting) initial state and then, after a sudden 
change of the interaction parameter U, evolves in time as pre¬ 
scribed by the interacting Hamiltonian. While the setup is 
apparently simple, the search for universal properties of the 
time evolution remains notoriously difficult due to the non- 
integrability of the Hubbard model in two and higher dimen¬ 
sions. Apart from the general assumption that non-integrable 
models feature thermalization and thus lose memory of the 
initial state in the long-time limit,!^ only the time evolution 
after quenches to a weak, finite Hubbard U seems to be well 
understood so far. Here, it coul d be show n by means of 
weak-coupling perturbation theotyEElIinEI] observables 
initially relax to non-thermal, quasistationary values (the sys¬ 
tem prethermalizes) before the significantly slower relaxation 
towards the thermal values sets in. 

It was later worked ou^ that the mechanism which traps 
the system in a quasi-stationary prethermal state is quite sim¬ 
ilar to the mechanism that hinders non-interacting systems 
from thermalizing. In the latter case the integrability of the 
Hamiltonian leads to a large number of constants of motion 
that highly constrain the dynamics of the system. In case 
of weakly interacting systems it is the proximity to the inte- 
grable point that introduces approximate constants of motion 
and hinders relaxation beyond the prethermalization plateau 
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on short timescales t <T/U'^ (here, T is the nearest-neighbor 
hopping). Relaxation towards the thermal average is delayed 
until later times (f > /U^). 

As a proof of concept of our formalism we use nonequilib¬ 
rium CPT to investigate the short- and long-time dynamics of 
an inhomogeneous initial state after an interaction quench in 
the Hubbard model. In particular we will study if and to what 
extent the CPT is able to describe prethermalization and the 
subsequent relaxation to a thermal state. 


B. Setup 

We consider the Hubbard model at zero temperature 
(j3 —> °°) and half-filling {{i = U /2) on a square lattice of 
L= 10 X 10 sites with periodic boundary conditions. Clus¬ 
ter indices run over 7,7 G {0,1,... ,24} and i,J G {0,1,2,3}, 
so that the system is cut into 25 clusters of size 2x2. The 
hopping is restricted to nearest neighbors and we set T = 1 
to fix energy and time units. Translational invariance of the 
initial state is broken by applying a local magnetic field of 
strength B to an arbitrarily chosen “impurity site” (here, site 0 
in cluster 0): 

T'/ait) = -ZaSij5ij5jfi5iflB{t) , (44) 

where ) is non-zero and unity for nearest neighbors only 
and where = 4-1 and Z| = — 1. Initially, the magnetic field 
is switched on with strength B{Q) = 10 to induce a (nearly) 
fully polarized magnetic moment on the impurity site and then 
switched off for times f > 0: 

B(f)=B(0)(l-©(f)). (45) 

Here, 0(f) is the Heaviside step function. Furthermore, the 
interaction t/(f) is switched off initially and then switched on 
to a non-zero value t/fin 

t/(f) = t/fin©(f). (46) 

Hence, in the quantum quench considered here, two parame¬ 
ters are changed simultaneously. The initial Hamiltonian 77ini 
features no interactions but is inhomogeneous due to the local 
magnetic field, the final Hamiltonian is translationally in¬ 
variant due to the absence of the magnetic field but has a finite 
interaction t/fln > 0. 

To apply nonequilibrium CPT, we use exact diagonalization 
to solve the 25 independent cluster problems and to construct 
the Hamiltonian of the effective medium (for details on the 
numerical implementation see Appendix [A|. Finally, Eq. ( |40l i 
is used to account for the inter-cluster hopping. The num¬ 
ber of non-zero elements of a cluster’s Q-matrix and therefore 
the computational effort of our approach increases quadrati- 
cally with the number of active states in the density matrix 
Pcluster ~ (^clusteri.e., 

states that contribute with a significant weight exp(—j3£m) to 
thermal averages. For convenience we have therefore chosen a 
zero-temperature initial state and consider a weak interaction 
U = 10^^ to lift the ground-state degeneracy present in the 



FIG. 4: (Color online) Time evolution of the local magnetic moment 
at the impurity (miinp(f), blue line) and its nearest neighbors (mNN(f), 
green line). The dark-blue (dark-green) arrow, which is pointing 
from right to left, indicates the long time average of the blue (green) 
curve. The light-blue (light-green) arrow, which is pointing from left 
to right, indicates the analytical average (|47}. The long time average 
was taken over 500,000 data points in the interval [0.5 x 10^, 10’*]. 


non-interacting system (denoted as U — 0^ in the following). 
The effective Hamiltonian h\t) for each cluster is then of size 
48 X 48 and the final CPT Hamiltonian of size 1200 x 1200. 
Exploiting its sparse form we are able to perform 1,000,000 
time steps with At = 0.01 to reach a maximal time fmax = 10^ 
with modest computational effort. Eor comparison we note 
that prior studies based on the nonequilibrium CPT, e.g. Refs. 
I7I22I have been limited to fmax = 10-20 inverse hoppings. 

The partitioning of the lattice into 2x2 clusters by CPT 
breaks rotational and reflection symmetries of the original 
problem. These are restored by averaging the resulting one- 
particle density matrix over the 4 possible ways to cut the 
lattice into 2x2 clusters. In the following we will show 
results for the time evolution of the local magnetic moment 
m,(f) = — n;|(f) at the impurity (mimp(f)) at its near¬ 

est neighbors (/wnn}?))- Only the latter are affected by the av¬ 
eraging. It restores the equivalence of nearest neighbors that 
lie in the same and nearest neighbors that lie in a neighboring 
cluster of the impurity. The extensive quantities total energy 
£’tot(0 = Ekin}?) +£’mt(t) (cf. Eq. ( |4^ and preceding discus¬ 
sion) and total magnetization M{t) = £, /«,(?) are both unaf¬ 
fected by the averaging. 

The initial state is the same for all quenches dis¬ 
cussed in the following. We And a polarization of 
'«imp(0) « 0.97 at the impurity which is partially screened 
(e.g., otnn(O) = —0.04) so that the total magnetization 
amounts to M(0) = « 0.70. 


C. Noninteracting case 

We first discuss the non-interacting case, i.e., a purely mag¬ 
netic quench where t/fin = 0+. Here, CPT predicts the ex¬ 
act time evolution (cf. the discussion below Eq. ( [37] l) since 
the cluster self-energies TI vanish. Our results are shown in 
Eig. 1^ Eor short times (t G [10^^,4 x 10**]) the local mag¬ 
netic moment at the impurity minip(f) (blue line) decays to a 
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t t t 

FIG. 5: (Color online) CPT results for the time evolution of the local magnetic moment at the impurity (minip(f), blue line) and at its nearest 
neighbors green line) for quenches from the limit of vanishing interaction (7 = 0+ (numerically implemented by setting U = 0.0001) 

to finite t/fin. In the insets the long-time behavior (f G [5 x 10^, 10^]) is plotted on a linear scale. The interval consists of 500,000 data points 
and was also used to calculate the long-time average (straight dashed lines). In total 1,000,000 time steps were performed with At = 0.01 on a 
L= 10 X 10 lattice (cut into 25 clusters of size 2 x 2 by CPT). 


value slightly above zero. Subsequently (f C [4 x 10°, 10^]) 
the dynamics is governed by collapse-and-revival oscillations 
caused by the finite system size. In particular we find that 
returns arbitrarily close to its initial value for large 
times. This is readily understood from the fact that the sys¬ 
tem’s dynamics is governed by the one-particle propagator 
exp(—iTfinO where Tf^^ denotes the final hopping matrix (i.e., 
after the quench). Tfln involves only a small number of differ¬ 
ent one-particle energy levels and thus U{t,Q) returns arbitrar¬ 
ily close to the identity matrix over time. 

For the non-interacting system it is possible to directly ac¬ 
cess the long-time average of the one-particle density matrix. 
One finds 

Pi;a= lim - - / dt{c]^{t)cja{t)) 

fmax JO 

kk' 


where we used that //fln can be diagonalized by a Fourier 
transformation involving the reciprocal lattice vectors k 
{Ri denotes the lattice vector to site i). We then have 

= I.la^Ajkc “d Cia{t) = 

where L is the system size. In Fig. ffl this prediction is com¬ 
pared with the numerical time average and indeed shows per¬ 
fect agreement. It is interesting to note that for non-degenerate 
energy levels one would have = Na/L, where Na 


is the total number of particles with spin a, and therefore 
=M{Q)/L. We conclude that degeneracy of energy lev¬ 
els is required to find memory of the initial state encoded in 
the average local magnetic moments mf'^. 


D. Quenches to finite (7fin 

For finite Upn CPT becomes an approximation and it is a 
priori unclear what kind of phenomena it is able to describe. 
In Fig. I^we show the long-time evolution for quenches to dif¬ 
ferent t/fin. For weak t/fin ^ 0.5 we find a (prethermalization- 
like) separation into two different time scales. Initially the 
time evolution qualitatively follows the non-interacting case, 
i.e., we see a fast decay of the local moment at the impu¬ 
rity site (blue line) followed by a quasi-stationary region of 
collapse-and-revival oscillations. For larger times these os¬ 
cillations decay and the system relaxes into a state character¬ 
ized by quasi-periodic fluctuations around its long-time aver¬ 
age (dashed blue line) which are driven by different frequen¬ 
cies. Taking a look at the t/fin dependence of the dynamics 
we notice that the region of collapse-and-revival oscillations 
shrinks with increasing t/fin and finally vanishes for t/fin ^ 1. 
The system then directly relaxes into a state with fluctuations 
around its long-time average. 

For comparison, also the magnetic moment at the neigh¬ 
bouring sites mNN(0 is plotted. While its dynamics for short 
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times must naturally be different from mimp (f) due to the inho¬ 
mogeneous initial state, we would expect a qualitative agree¬ 
ment in the long-time limit if the system thermalizes. How¬ 
ever, this is not the case. There remains a clear difference in 
the amplitude of the fluctuations around the long-time average 
up to the largest simulated times. Hence we conclude that the 
system still keeps memory of the initial state and thus does 
not thermalize. 

Having in mind the general discussion on prethermaliza- 
tion in Sec. V A one can give an intuitive interpretation of 
these observations based on the effective-medium approach: 
While the non-interacting system is isolated and its dynamics 
is constrained through many constants of motion, there is a 
large number of virtual orbitals coupled to the system in the 
interacting case. These virtual orbitals act like a surround¬ 
ing bath. For weak t/gn the virtual orbitals are only weakly 
coupled to the system and their influence is delayed to large 
times, while initially the dynamics is constrained similar to 
the non-interacting case. For strong f/fin, on the other hand, 
the coupling is strong and affects the dynamics of the sys¬ 
tem considerably. However, the number of virtual sites is still 
too small to allow for a complete dissipation of the informa¬ 
tion on the initial state into the bath. Therefore, a thermalized 
state is not reached. For an exact calculation the number of 
virtual sites would scale exponentially in system size. For 
CPT, on the other hand, it scales exponentially only in clus¬ 
ter size but linearly in the number of clusters and thus in the 
system size. Memory of the initial state is therefore retained 
within the one-particle density matrix and leaves its traces in 
the magnetic moments as seen in our calculations. 


E. Violation of conservation laws 



-280 . Vtn . . .. I I . . r . I m ill 

10-2 jQ-l igO jQl jq2 jq3 jo4 


t 

FIG. 6: (Color online) Violation of energy conservation by CPT. The 
numbers indicate the respective value of (/(;„. Energy conservation is 
respected for t/fin = 0+ where CPT is exact (blue line). An increas¬ 
ingly significant violation of energy conservation is seen for larger 
Ffin. 



FIG. 7: (Color online) Violation of conservation of total magnetiza¬ 
tion M by CPT. The numbers indicate the value of (/flu. Curves for 
Ffin > 0.25 are only partially plotted for better visibility. 


CPT as an approximation lacks any kind of self-consistency 
and is thus unable to respect the fundamental continuity equa¬ 
tions and their corresponding conservation laws.f^ Therefore, 
one has to expect a violation of energy- or particle-number 
conservation, for example. Furthermore, in contrast to the 
equilibrium case where CPT interpolates between the ex¬ 
act limits U = 0 and T = 0, it yields exact results only for 
quenches to Ufm = 0. The dynamics after a quench to the 
atomic limit Tfin = 0 (with finite Unn > 0) cannot be described 
exactly due to the non-local entanglement of the initial state. 
We thus generally expect that the quality of the CPT results 
degrades with increasing interaction strength. 

The numerical results for the total energy, see Fig. con¬ 
firm this expectation. Energy conservation is respected for 
Ufin = 0, where CPT is exact. With t/fin > 0 and increasing, 
however, a significant time dependence of the total energy sets 
in earlier and earlier. For t/fin ^ 1 energy conservation is vio¬ 
lated already for f < 10. Similar results are found for the total 
magnetization M — Li(^I't ~”ii)’ <2f. Fig. While the magne¬ 
tization should be constant for all times since neither hopping 
nor interaction (cf. Eqs. ( |44l l and ( |46l l) involve spin-flip terms, 
we And such behavior only for short times. Eor longer times 
oscillations arise and the conservation of total magnetization 
is violated. Eor increasing t/fin the oscillations set in earlier 


indicating again that the quality of CPT is best for values of 
t/fin close to zero. 

We note that the total particle number A = -F A|, how¬ 
ever, is conserved during the time evolution. This holds true 
for a half-filled and homogeneously charged system and is due 
to the fact that CPT preserves particle-hole symmetry. This 
can easily be understood as follows: Each cluster Hamiltonian 
is particle-hole symmetric and since each cluster is solved 
exactly within CPT the corresponding effective Hamiltonian 
h\t) is also particle-hole symmetric. The CPT Hamiltonian 
is now given by Eq. ( |40l l which additionally includes the inter¬ 
cluster hopping. However, the inter-cluster hopping is clearly 
particle-hole symmetric and so is the final CPT Hamiltonian. 

VI. SUMMARY AND OUTLOOK 

Concluding, we have shown that the nonequilibrium self¬ 
energy of an interacting lattice-fermion model can uniquely 
be decomposed into a superposition of noninteracting, iso¬ 
lated modes. This decomposition is a direct analog to a well- 
established decomposition of equilibrium Green’s functions, 
called the Lehmann representation. Our proof not only pro¬ 
vides a direct scheme to construct the Lehmann representation 
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of the self-energy, and thus allows for a deeper theoretical un¬ 
derstanding of the self-energy complementary to its diagram¬ 
matic definition, but also proves useful for practical applica¬ 
tions. 

As a proof of concept we investigated the time evolution 
of local magnetic moments in the fermionic Hubbard model 
after an interaction quench using nonequilibrium cluster- 
perturbation theory. Our formalism allowed to avoid the so¬ 
lution of an inhomogeneous Dyson equation on the Keldysh 
contour and we were able to propagate the one-particle den¬ 
sity matrix up to times fmax = lO'*. 

On the physical side, quenches to weak Ufm turned out to be 
most interesting. In agreer nent with the predictions of general 
perturbative considerations,ElllE3E2|yYe found a separation of 
the dynamics into two time scales. While the system qualita¬ 
tively follows the constrained dynamics of the non-interacting 
t/fln = 0 limit, the constraints are broken up for large times due 
to the interaction and the system shows signs of relaxation. 
However, memory of the initial state persists in the density 
matrix up to the largest simulated times clearly indicating the 
absence of thermalization. 

While the simple treatment of correlations by nonequilib¬ 
rium CPT has shown to be enough to cover the mentioned 
two-stage relaxation dynamics, it also leads to a violation of 
the fundamental conservation laws of energy and total magne¬ 
tization. This could be fixed by additionally imposing a self- 
consistency condition as it is done in nonequilibrium DMFT 
or in self-energy functional theory. Due to the significant, ad¬ 
ditional complexity of these approaches, however, simulations 
would again be restricted to short time scales. A simpler, more 
pragmatic approach might thus be preferable where, for exam¬ 
ple, local continuity equations are enforced to ensure energy, 
total magnetization and particle-number conservation.QSI Such 
a “conserving cluster-perturbation theory” could allow for a 
complete dissipation of initial perturbations and thus total loss 
of the memory of the initial state. Work along these lines is in 
progress. 
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Appendix A: Numerical construction of the effective 
Hamiltonian 

1. The g-matrix and its time derivatives 

We assume that a small cluster is solved using exact diag- 
onalization and that all time derivatives = d"H{t) of 

the Hamiltonian are known analytically. The numerical eval¬ 
uation of Eq. Q for the g-matrix is straightforward within 
exact diagonalization. Its n-th derivative can be obtained as 
follows. We have 

t/W(f,o) = (Ai) 

k=o\ ^ / 

for the propagator t/(f,0). The n-th derivative can 

then be calculated iteratively as it only depends on 
with k < n. Using further that 



one finds the n-th derivative cf‘\t) of the annihilation opera¬ 
tor and thus of see Eq. 0- In the following we will 

assume that is available to arbitrary order. 

2. Construction of the effective Hamiltonian at f = 0 

We start by constructing 2^(0), i.e., a basis for the virtual 
sector. It is easy to verify that 

Paa'=Y.Q*a{mia'{0), (A3) 

i 

defines a projector. Diagonalization of P yields the eigenval¬ 
ues 0 and 1. Eigenvectors corresponding to 1 are given by 
2(0)^ itself, eigenvectors corresponding to 0 form the desired 
matrix [2^(0)]^. Initially, the effective medium is in equilib¬ 
rium and thus explicitly given by Eq. ( [TtI i at f = 0. However, 
since we picked the completing basis vectors arbitrarily, we 
will have ^ 0 for s ^ s\ i.e., generally h will not be diag¬ 
onal in the virtual sector. Explicit diagonalization of h in the 
virtual sector yields a unitary transform R 

K,,=Y^RsrdrR;,,. (A4) 

r 

Replacing 2^(0) ^2^(0)^ we get —>■ 5^^ids, i.e., we 

have found a completing basis so that h is diagonal in the vir¬ 
tual sector. 

3. The time derivatives (f) 

Assume that h{t),Q{t),Q^{t) and ^re known for 

an arbitrary time t. This is at least the case for f = 0 as we have 
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seen so far. We recall that we required !i{t) to be constant in 0{t +At). Analytically this can be written as 
the virtual sector (cf. discussion below Eq. ([H) 


/Ay(0 = 5,.,'/r..(0) ^ (A5) 

i.e., all time derivatives vanish in the virtual sector. Only 
the hybridization elements and the physical sector yield non¬ 
trivial elements. They follow from Eq. as 

t (l) 

k=o a 

(A6) 

(f) on the other hand only depends on (t), and (f), 
for k <n, as readily follows from 

= (A7) 

= -'L 

It is thus possible to iteratively calculate (9^") (f) and (f). 


0(f+Af)=r 



/ rt+At 

{-•I ' 


h{t')dt'] \o{t). 


(A8) 


Using the Magnus expansion EH the propagator can be system¬ 
atically expanded in Af” and Assuming that Af lies 

within the convergence radius of the Magnus expansion (this 
is generally expected to be the case for sufficiently small Af), 
we can reduce the propagation error arbitrarily by increasing 
the order. In practice, an evaluation of the Magnus expan¬ 
sion using commutator-free exponential time propagatorJSl 
(CEETs) allows for an efficient numerical propagation which 
takes advantage of the sparse form of the effective Hamilto¬ 
nian. 

Having found (9(f -f Af), we get h{t + At) from 


/f,y(f+Af) =fXeW(f+Af)e-‘^“'[0(f + Af)1']„,, (A9) 

a 


4. Propagation of the 0-matrix 

We assume that 0(f) and all derivatives of (f) are known 
at some time f and we want to propagate the 0-matrix to 


and can thus proceed by calculating O^”^ (f + At) and 
/;(”)(f _l_ Af) completing the circle. We emphasize that the 
whole procedure is numerically exact, i.e., the error is below 
machine precision, if Af is chosen sufficiently small. 
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